library(mvtnorm)
library(tmvtnorm)
library(bayesm)
dyn.load("test/Orthant.so")

A<-matrix(c(1,-0.95,-0.97,-0.94,0.97,-0.95,1,-0.99,-0.99,-0.99,0.97,0.99,1,-0.99,0.99,-0.94,-0.99,-0.99,1,-0.99,-0.97,-0.99,0.99,-0.99,1),nrow=5)
a<-rep(0.05,p)

ft<-function(S,a,m,M)
{
	p<-ncol(S)
	theta<-0
	res1<-.C("sampler", as.integer(M),as.integer(p),as.double(a), as.double(S),as.double(m),res=as.double(theta))
	return(res1$res)
}

it<-function(S,a,m,M)
{
	p<-ncol(S)
	theta<-0
	res1<-.C("sampler2", as.integer(M),as.integer(p),as.double(a), as.double(S),as.double(m),res=as.double(theta))
	return(res1$res)
}
kt<-function(S,a,m,M)
{
	p<-ncol(S)
	L<-t(chol(S))
	theta<-0
	res1<-.C("ghk_oneside", as.double(L),as.double(a-m),as.double(rep(0,p)),as.integer(p),as.integer(M),res=as.double(theta))
	return(res1$res)
}


Sub<-function(A, i ,j)
{
	p<-nrow(as.matrix(A))	
	if(i>j)
	{
		s<-j
		j<-i
		i<-s
	}
	Y<-matrix(.C("Substitute",as.double(A),res=as.double(A),as.integer(p),as.integer(i),as.integer(j))$res,nrow=p)
}
order<-function(S,a,m,M)
{
	pp<-ncol(as.matrix(S))
	a<-a-m
	for(i in 1:pp)
	{	
		resp<-1
		for( j in i:pp)
		{
			SS<-Sub(S,i,j)
			SS<-SS[1:i,1:i]
			aa<-a
			aa[i]<-a[j]
			aa<-aa[1:i]
			p<-ncol(as.matrix(SS))
			L<-t(chol(SS))
			theta<-0
			res<-.C("ghk_order", as.double(L),as.double(aa),as.double(rep(0,p)),as.integer(p),res=as.double(theta))$res
			if(res<resp)
			{
				resp<-res
				k<-j
			}		
		}
			S<-Sub(S,i,k)
			s<-a[i]
			a[i]<-a[j]
			a[j]<-s
	}
	return(list(S=S,a=a))
}

jt<-function(S,a,m,M,tr)
{
	p<-ncol(S)
	L<-t(chol(S))
	theta<-0
	res1<-.C("sampler3", as.integer(M),as.integer(p), as.double(L),as.double(a-m),as.double(rep(0,p)),as.double(0.9), as.integer(tr),res=as.double(theta))
	return(res1$res)
}

ist<-function(S,a,m,M)
{
	p<-ncol(S)
	L<-t(chol(S))
	theta<-0
	res1<-.C("IS", as.integer(M),as.integer(p), as.double(L),as.double(a-m),res=as.double(theta))
	return(res1$res)
}
mt<-function(S,a,m,M,d)
{
	p<-ncol(S)
	L<-t(chol(S))
	theta<-0
	res1<-.C("sampler4", as.integer(M),as.integer(p), as.double(L),as.double(a-m),as.double(rep(0,p)),as.double(0.9),as.intger(d),res=as.double(theta))
	return(res1$res)
}

mst<-function(S,a,m,M)
{
	p<-ncol(S)
	L<-t(chol(S))
	theta<-0
	nu<-3
	res1<-.C("sampler5", as.integer(M),as.integer(p), as.double(L),as.double(a-m),as.integer(nu),as.double(0.9),res=as.double(theta))
	return(res1$res)
}
gt<-function(S,a,m)
{
	res<-pmvnorm(lower=a,upper=Inf,mean=m,sigma=S)
	return(res)
}

ht<-function(S,a,m,M)
{
	C<-t(chol(S))
	p<-length(a)
	res<-ghkvec(L=C,trunpt=a-m,above=rep(0,p),r=M)
	return(res)
}


